/*
 * Test.cpp
 *
 *  Created on: 23 May 2011
 *      Author: allan
 */

// C++ includes
#include <cmath>
#include <iostream>
#include <iomanip>
using namespace std;

// FluidReact includes
#include "AqueousPhase.h"
#include "System.h"
#include "Speciation.h"
#include "Reaction.h"
#include "BatchReactor.h"

void TestActivityCoefficients_NaCl()
{
	AqueousPhase aqueousPhase;
	
	aqueousPhase.AddSpecies("H2O(l)");
	aqueousPhase.AddSpecies("H+");
	aqueousPhase.AddSpecies("OH-");
	aqueousPhase.AddSpecies("Na+");
	aqueousPhase.AddSpecies("Cl-");
	aqueousPhase.AddSpecies("NaCl(aq)");
	
	aqueousPhase.SetElectrolyte("NaCl(aq)", vector<pair<string,double>>({{"Na+", +1.0}, {"Cl-", -1.0}}));
	
	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
	
	cout << setw(15) << "MolalityNaCl";
	cout << setw(15) << "IonicStrength";
	cout << setw(15) << "ActivityH2O";
	cout << setw(15) << "GammaNaCl+-";
	cout << setw(15) << "GammaH[+]";
	cout << setw(15) << "GammaOH[-]";
	cout << setw(15) << "GammaNa[+]";
	cout << setw(15) << "GammaCl[-]";
	cout << setw(15) << "GammaNaCl(aq)" << endl;
		
	for(double mNaCl = 0.1; mNaCl < 6.0; mNaCl += 0.1)
	{
		aqueousPhase.SetMolesOfSpecies({55.508, 0.0, 0.0, mNaCl, mNaCl, 0.0});
		
		vector<double> activityCoeffs = aqueousPhase.ActivityCoefficients(25.0, 100.0);
		
		cout << fixed << showpoint;
		
		cout << setw(15) << mNaCl;
		cout << setw(15) << aqueousPhase.GetEffectiveIonicStrength();
		cout << setw(15) << activityCoeffs[0] * 55.508;
		cout << setw(15) << sqrt(activityCoeffs[3] * activityCoeffs[4]);
		cout << setw(15) << activityCoeffs[1];
		cout << setw(15) << activityCoeffs[2];
		cout << setw(15) << activityCoeffs[3];
		cout << setw(15) << activityCoeffs[4];
		cout << setw(15) << activityCoeffs[5];
		cout << endl;
	}
}

void TestSpeciationWaterNaCl()
{
	AqueousPhase aqueousPhase;
	
	aqueousPhase.AddSpecies("H2O(l)");
	aqueousPhase.AddSpecies("H+");
	aqueousPhase.AddSpecies("OH-");
	aqueousPhase.AddSpecies("Na+");
	aqueousPhase.AddSpecies("Cl-");
	aqueousPhase.AddSpecies("NaCl(aq)");
	
	aqueousPhase.SetElectrolyte("NaCl(aq)", vector<pair<string,double>>({{"Na+", +1.0}, {"Cl-", -1.0}}));
	
	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);

	System system;

	system.SetTemperature(25.0);
	system.SetPressure(100.0);
	
	system.AddPhase(aqueousPhase);

	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-"};
	vector<string> equilibriumSpecies = {"NaCl(aq)", "OH-"};
	
	system.SetPrimarySpecies(primarySpecies);
	system.SetEquilibriumSpecies(equilibriumSpecies);
	
	EquilibriumReaction eqReactionNaClaq;
	
	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
	
	EquilibriumReaction eqReactionOH;
	
	eqReactionOH.SetNameOfReactant("OH-");
	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
	
	system.AddEquilibriumReaction(eqReactionNaClaq);
	system.AddEquilibriumReaction(eqReactionOH);
	
//	double C_H2O  = 55508;
//	double C_H    = 1.0E-7;
//	double C_OH   = 1.0E-7;
//	double C_Na   = 1.0E-7;
//	double C_Cl   = 1.0E-7;
//	double C_NaCl = 1000;
	
	double C_H2O  = 55508;
	double C_H    = 1.0E-4;
	double C_OH   = 1.0E-4;
	double C_Na   = 960;
	double C_Cl   = 960;
	double C_NaCl = 40;
	
	system.SetMolesOfPhase(0, {C_H2O, C_H, C_OH, C_Na, C_Cl, C_NaCl});
	
	Speciation speciation;
	
	speciation.SetSystem(system);
	
	speciation.ImposeUj(0, 55.508);
	speciation.ImposeUj(1, 0.0);
	speciation.ImposeUj(2, 1000);
	//speciation.ImposeUj(3, U_Cl);
	
	speciation.ImposeCe(0, 39);
	
	speciation.Solve(1.0E-6, 100);
	
	cout << system.GetMolesOfSpecies(0) << endl;
	cout << system.GetMolesOfSpecies(1) << endl;
	cout << system.GetMolesOfSpecies(2) << endl;
	cout << system.GetMolesOfSpecies(3) << endl;
	cout << system.GetMolesOfSpecies(4) << endl;
	cout << system.GetMolesOfSpecies(5) << endl;
	
	cout << endl;
	
	cout << system.GetActivityOfSpecies(0) << endl;
	cout << system.GetActivityOfSpecies(1) << endl;
	cout << system.GetActivityOfSpecies(2) << endl;
	cout << system.GetActivityOfSpecies(3) << endl;
	cout << system.GetActivityOfSpecies(4) << endl;
	cout << system.GetActivityOfSpecies(5) << endl;
}

void TestBatchReactorWaterNaClCaCO3()
{
	AqueousPhase aqueousPhase;
	
	aqueousPhase.AddSpecies("H2O(l)");
	aqueousPhase.AddSpecies("H+");
	aqueousPhase.AddSpecies("OH-");
	aqueousPhase.AddSpecies("Na+");
	aqueousPhase.AddSpecies("Cl-");
	aqueousPhase.AddSpecies("Ca++");
	aqueousPhase.AddSpecies("HCO3-");
	aqueousPhase.AddSpecies("NaCl(aq)");
	aqueousPhase.AddSpecies("CaCO3(aq)");
	
	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
	aqueousPhase.SetActivityCoefficientModelSetschenow("CaCO3(aq)", 0.10);

	MineralPhase calcitePhase;
	
	calcitePhase.AddSpecies("Calcite");
	
	EquilibriumReaction eqReactionOH;
	
	eqReactionOH.SetNameOfReactant("OH-");
	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
	
	EquilibriumReaction eqReactionNaClaq;
	
	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
	
	EquilibriumReaction eqReactionCaCO3aq;
	
	eqReactionCaCO3aq.SetNameOfReactant("CaCO3(aq)");
	eqReactionCaCO3aq.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
	eqReactionCaCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	eqReactionCaCO3aq.SetEquilibriumConstant({7.09941}, {25.0});
	
	MineralReaction minReactionCalcite;
	
	minReactionCalcite.SetNameOfReactant("Calcite");
	minReactionCalcite.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
	minReactionCalcite.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	minReactionCalcite.SetEquilibriumConstant({1.84377}, {25.0});
	minReactionCalcite.SetSpecificSurfaceArea(9.8 * 100.09); // 9.8 cm2/g * 100.09 g/mol
	
	MineralReactionMechanism calciteMechanism;
	
	calciteMechanism.k  = pow(10.0, -3.48);
	calciteMechanism.Ea = 35.4;
	
	minReactionCalcite.AddMechanisms(calciteMechanism);

	System system;

	system.AddPhase(aqueousPhase);
	system.AddPhase(calcitePhase);

	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "Ca++", "HCO3-"};
	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "CaCO3(aq)"};
	vector<string> kineticSpecies = {"Calcite"};
	
	system.SetPrimarySpecies(primarySpecies);
	system.SetEquilibriumSpecies(equilibriumSpecies);
	system.SetKineticSpecies(kineticSpecies);
	
	system.AddEquilibriumReaction(eqReactionOH);
	system.AddEquilibriumReaction(eqReactionNaClaq);
	system.AddEquilibriumReaction(eqReactionCaCO3aq);
	
	system.AddKineticReaction(minReactionCalcite);
	
	double C_H2O      = 55508;
	double C_H        = 1.0E-7;
	double C_OH       = 1.0E-4;
	double C_Na       = 1.0E-4;
	double C_Cl       = 1.0E-4;
	double C_Ca       = 1.0E-4;
	double C_HCO3     = 1.0E-4;
	double C_NaClaq   = 1000;
	double C_CaCO3aq  = 1.0E-4;
	
	double C_Calcite  = 100;
	
	system.SetMolesOfPhase(0, {C_H2O, C_H, C_OH, C_Na, C_Cl, C_Ca, C_HCO3, C_NaClaq, C_CaCO3aq});
	system.SetMolesOfPhase(1, {C_Calcite});
	
	BatchReactor batchReactor;
	
	batchReactor.SetSystem(&system);
	
	batchReactor.Solve(3, 1.0E-5, 1.0E-3, "BatchReactor[Water-NaCl-Calcite][Simulation 1].dat");
}

void TestSpeciationWaterCO2()
{
	AqueousPhase aqueousPhase;
	
	aqueousPhase.AddSpecies("H2O(l)");
	aqueousPhase.AddSpecies( "H+");
	aqueousPhase.AddSpecies("OH-");
	aqueousPhase.AddSpecies("HCO3-");
	aqueousPhase.AddSpecies("CO2(aq)");
	
	GaseousPhase gaseousPhase;
	
	gaseousPhase.AddSpecies("H2O(g)");
	gaseousPhase.AddSpecies("CO2(g)");
	
	EquilibriumReaction eqReactionOH;
	
	eqReactionOH.SetNameOfReactant("OH-");
	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
	
	EquilibriumReaction eqReactionCO2aq;
	
	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2aq.SetEquilibriumConstant({-6.35332}, {25.0});
	
	EquilibriumReaction eqReactionH2Og;
	
	eqReactionH2Og.SetNameOfReactant("H2O(g)");
	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
	eqReactionH2Og.SetEquilibriumConstant({1.49929}, {25.0});
	
	EquilibriumReaction eqReactionCO2g;
	
	eqReactionCO2g.SetNameOfReactant("CO2(g)");
	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2g.SetEquilibriumConstant({-7.82545}, {25.0});
	
	System system;

	system.SetTemperature(25.0);
	system.SetPressure(100.0);
	
	system.AddPhase(aqueousPhase);
	system.AddPhase(gaseousPhase);

	vector<string> primarySpecies = {"H2O(l)", "H+", "HCO3-"};
	vector<string> equilibriumSpecies = {"OH-", "CO2(aq)", "H2O(g)", "CO2(g)"};
	
	system.SetPrimarySpecies(primarySpecies);
	system.SetEquilibriumSpecies(equilibriumSpecies);
	
	system.AddEquilibriumReaction(eqReactionOH);
	system.AddEquilibriumReaction(eqReactionCO2aq);
	system.AddEquilibriumReaction(eqReactionH2Og);
	system.AddEquilibriumReaction(eqReactionCO2g);
	
	double C_H2O      = 55508;
	double C_H        = 1.0E-7;
	double C_OH       = 1.0E-4;
	double C_HCO3     = 1.0E-4;
	double C_CO2aq    = 1.0E-4;
	
	double C_H2Og     = 1.0E-4;
	double C_CO2g     = 5;
	
	system.SetMolesOfPhase(0, {C_H2O, C_H, C_OH, C_HCO3, C_CO2aq});
	system.SetMolesOfPhase(1, {C_H2Og, C_CO2g});
	
	Speciation speciation;
	
	speciation.SetSystem(system);
	
	speciation.ImposeChargeBalance();
	speciation.ImposeCj(0, 55.508);
	speciation.ImposeCe(3, 30);
	
	speciation.Solve(1.0E-6, 100);
	
	cout << system.GetMolesOfSpecies(0) << endl;
	cout << system.GetMolesOfSpecies(1) << endl;
	cout << system.GetMolesOfSpecies(2) << endl;
	cout << system.GetMolesOfSpecies(3) << endl;
	cout << system.GetMolesOfSpecies(4) << endl;
	cout << system.GetMolesOfSpecies(5) << endl;
	cout << system.GetMolesOfSpecies(6) << endl;
	
	cout << endl;
	
	cout << system.GetActivityOfSpecies(0) << endl;
	cout << system.GetActivityOfSpecies(1) << endl;
	cout << system.GetActivityOfSpecies(2) << endl;
	cout << system.GetActivityOfSpecies(3) << endl;
	cout << system.GetActivityOfSpecies(4) << endl;
	cout << system.GetActivityOfSpecies(5) << endl;
	cout << system.GetActivityOfSpecies(6) << endl;
	
	cout << endl;
	
	cout << log10(system.GetEquilibriumReactionQuotient(0)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(1)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(2)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(3)) << endl;
	
	cout << endl;
	
	cout << "I  = " << system.GetAqueousPhase()->GetEffectiveIonicStrength() << endl;
	cout << "aw = " << system.GetAqueousPhase()->Activities(25.0, 100.0)[0] << endl;
	
	cout << endl;
	
	cout << system.GetMolesOfSpecies(1) - system.GetMolesOfSpecies(2) - system.GetMolesOfSpecies(3) << endl;
}

void TestSpeciationWaterNaClCO2()
{
	AqueousPhase aqueousPhase;
	aqueousPhase.AddSpecies("H2O(l)");
	aqueousPhase.AddSpecies("H+");
	aqueousPhase.AddSpecies("OH-");
	aqueousPhase.AddSpecies("Na+");
	aqueousPhase.AddSpecies("Cl-");
	aqueousPhase.AddSpecies("HCO3-");
	aqueousPhase.AddSpecies("NaCl(aq)");
	aqueousPhase.AddSpecies("HCl(aq)");
	aqueousPhase.AddSpecies("NaOH(aq)");
	aqueousPhase.AddSpecies("NaHCO3(aq)");
	aqueousPhase.AddSpecies("CO2(aq)");
	aqueousPhase.AddSpecies("CO3--");

	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
	
	GaseousPhase gaseousPhase;
	gaseousPhase.AddSpecies("H2O(g)");
	gaseousPhase.AddSpecies("CO2(g)");
	
	EquilibriumReaction eqReactionOH;
	
	eqReactionOH.SetNameOfReactant("OH-");
	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
	eqReactionOH.SetEquilibriumConstant({14.0013, 12.03594}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionNaClaq;
	
	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionNaClaq.SetEquilibriumConstant({0.7769, 0.473}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionHClaq;
	
	eqReactionHClaq.SetNameOfReactant("HCl(aq)");
	eqReactionHClaq.SetNameOfProducts(vector<string>({"H+", "Cl-"}));
	eqReactionHClaq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionHClaq.SetEquilibriumConstant({0.69754, 0.62}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionNaOHaq;
	
	eqReactionNaOHaq.SetNameOfReactant("NaOH(aq)");
	eqReactionNaOHaq.SetNameOfProducts(vector<string>({"Na+", "H2O(l)", "H+"}));
	eqReactionNaOHaq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	eqReactionNaOHaq.SetEquilibriumConstant({14.18, 12.19733}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionNaHCO3aq;
	
	eqReactionNaHCO3aq.SetNameOfReactant("NaHCO3(aq)");
	eqReactionNaHCO3aq.SetNameOfProducts(vector<string>({"Na+", "HCO3-"}));
	eqReactionNaHCO3aq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionNaHCO3aq.SetEquilibriumConstant({0.24474, 0.72845}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionCO2aq;
	
	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2aq.SetEquilibriumConstant({-6.35332, -6.03078}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionCO3;
	
	eqReactionCO3.SetNameOfReactant("CO3--");
	eqReactionCO3.SetNameOfProducts(vector<string>({"HCO3-", "H+"}));
	eqReactionCO3.SetStoichiometryOfProducts({1.0, -1.0});
	eqReactionCO3.SetEquilibriumConstant({10.32666, 9.80904}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionH2Og;
	
	eqReactionH2Og.SetNameOfReactant("H2O(g)");
	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
	eqReactionH2Og.SetEquilibriumConstant({1.49929, -0.05018}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionCO2g;
	
	eqReactionCO2g.SetNameOfReactant("CO2(g)");
	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2g.SetEquilibriumConstant({-7.82545, -8.19835}, {25.0, 100.0});
	
	System system;

	system.SetTemperature(100.0);
	system.SetPressure(100.0);
	
	system.AddPhase(aqueousPhase);
	system.AddPhase(gaseousPhase);

	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "HCO3-"};
	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "HCl(aq)", "NaOH(aq)", "NaHCO3(aq)", "CO2(aq)", "CO3--", "H2O(g)", "CO2(g)"};
	vector<string> kineticSpecies = {};
	
	system.SetPrimarySpecies(primarySpecies);
	system.SetEquilibriumSpecies(equilibriumSpecies);
	system.SetKineticSpecies(kineticSpecies);
	
	system.AddEquilibriumReaction(eqReactionOH);
	system.AddEquilibriumReaction(eqReactionNaClaq);
	system.AddEquilibriumReaction(eqReactionHClaq);
	system.AddEquilibriumReaction(eqReactionNaOHaq);
	system.AddEquilibriumReaction(eqReactionNaHCO3aq);
	system.AddEquilibriumReaction(eqReactionCO2aq);
	system.AddEquilibriumReaction(eqReactionCO3);
	system.AddEquilibriumReaction(eqReactionH2Og);
	system.AddEquilibriumReaction(eqReactionCO2g);
	
	vector<double> C(system.GetNumOfSpecies(), 1.0E-7);
	
	C[system.GetIndexOfSpecies("H2O(l)")]    = 55.508;
	C[system.GetIndexOfSpecies("Na+")]       = 1.0;
	C[system.GetIndexOfSpecies("Cl-")]       = 1.0;
	C[system.GetIndexOfSpecies("HCO3-")]     = 2.0;
	C[system.GetIndexOfSpecies("CO2(g)")]    = 2.0; // Moleweight: 44.01 g/mol
	
	system.SetMolesOfSpecies(C);
	
	Speciation speciation;
	
	speciation.SetSystem(system);
	
	speciation.ImposeChargeBalance();
	speciation.ImposeCj(0, 55.508);	// H2O
	speciation.ImposeCj(2, 1.0);	// Na+
	speciation.ImposeCj(3, 1.0);	// Cl-
	speciation.ImposeCe(8, 2.0);	// CO2(g) (1 mole = 44 grams)
	
	speciation.Solve(1.0E-6, 1000);
	
	cout << endl;
	
	for(unsigned iSpecies = 0; iSpecies < system.GetNumOfSpecies(); ++iSpecies)
		cout << setw(10) << system.GetSpecies(iSpecies).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(iSpecies)  << endl;
	
	cout << endl;
	
	for(unsigned iSpecies = 0; iSpecies < system.GetNumOfSpecies(); ++iSpecies)
		cout << setw(10) << system.GetSpecies(iSpecies).GetName()  << ": " << setw(15) << system.GetActivityOfSpecies(iSpecies)  << endl;
	
	cout << endl;
	
	for(unsigned iEquilibrium = 0; iEquilibrium < system.GetNumOfEquilibriumSpecies(); ++iEquilibrium)
		cout << "K[" << system.GetSpecies(system.GetIndexOfEquilibriumSpecies(iEquilibrium)).GetName()  << "]: " << setw(15) << log10(system.GetEquilibriumReactionQuotient(iEquilibrium)) << endl;
}

void TestSpeciationWaterNaClCaCO3CO2()
{
	AqueousPhase aqueousPhase;
	
	aqueousPhase.AddSpecies("H2O(l)");
	aqueousPhase.AddSpecies("H+");
	aqueousPhase.AddSpecies("OH-");
	aqueousPhase.AddSpecies("Na+");
	aqueousPhase.AddSpecies("Cl-");
	aqueousPhase.AddSpecies("Ca++");
	aqueousPhase.AddSpecies("HCO3-");
	aqueousPhase.AddSpecies("NaCl(aq)");
	aqueousPhase.AddSpecies("CO2(aq)");
	aqueousPhase.AddSpecies("CaCO3(aq)");
	
	GaseousPhase gaseousPhase;
	
	gaseousPhase.AddSpecies("H2O(g)");
	gaseousPhase.AddSpecies("CO2(g)");
	
	EquilibriumReaction eqReactionOH;
	
	eqReactionOH.SetNameOfReactant("OH-");
	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
	
	EquilibriumReaction eqReactionNaClaq;
	
	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
	
	EquilibriumReaction eqReactionCO2aq;
	
	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2aq.SetEquilibriumConstant({-6.35332}, {25.0});
	
	EquilibriumReaction eqReactionCaCO3aq;
	
	eqReactionCaCO3aq.SetNameOfReactant("CaCO3(aq)");
	eqReactionCaCO3aq.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
	eqReactionCaCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	eqReactionCaCO3aq.SetEquilibriumConstant({7.09941}, {25.0});
	
	EquilibriumReaction eqReactionH2Og;
	
	eqReactionH2Og.SetNameOfReactant("H2O(g)");
	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
	eqReactionH2Og.SetEquilibriumConstant({1.49929}, {25.0});
	
	EquilibriumReaction eqReactionCO2g;
	
	eqReactionCO2g.SetNameOfReactant("CO2(g)");
	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2g.SetEquilibriumConstant({-7.82545}, {25.0});
	
	System system;

	system.SetTemperature(25.0);
	system.SetPressure(100.0);
	
	system.AddPhase(aqueousPhase);
	system.AddPhase(gaseousPhase);

	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "Ca++", "HCO3-"};
	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "CO2(aq)", "CaCO3(aq)", "H2O(g)", "CO2(g)"};
	
	system.SetPrimarySpecies(primarySpecies);
	system.SetEquilibriumSpecies(equilibriumSpecies);
	
	system.AddEquilibriumReaction(eqReactionOH);
	system.AddEquilibriumReaction(eqReactionNaClaq);
	system.AddEquilibriumReaction(eqReactionCO2aq);
	system.AddEquilibriumReaction(eqReactionCaCO3aq);
	system.AddEquilibriumReaction(eqReactionH2Og);
	system.AddEquilibriumReaction(eqReactionCO2g);
	
	double C_H2O      = 55.508;
	double C_H        = 1.0E-4;
	double C_OH       = 1.0E-4;
	double C_Na       = 1.0E-4;
	double C_Cl       = 1.0E-4;
	double C_Ca       = 1.0E-4;
	double C_HCO3     = 1.0E-4;
	double C_NaClaq   = 1.0E-4;
	double C_CO2aq    = 1.0E-4;
	double C_CaCO3aq  = 1.0E-4;
	
	double C_H2Og     = 1.0E-4;
	double C_CO2g     = 30;
	
	system.SetMolesOfPhase(0, {C_H2O, C_H, C_OH, C_Na, C_Cl, C_Ca, C_HCO3, C_NaClaq, C_CO2aq, C_CaCO3aq});
	system.SetMolesOfPhase(1, {C_H2Og, C_CO2g});
	
	Speciation speciation;
	
	speciation.SetSystem(system);
	
	speciation.ImposeChargeBalance();
	speciation.ImposeCj(0, 55.508);
	speciation.ImposeUj(2, 1);
	speciation.ImposeUj(3, 1);
	speciation.ImposeUj(4, 1);
	speciation.ImposeCe(5, 100);
	
	speciation.Solve(1.0E-6, 100);
	
	cout << setw(10) << system.GetSpecies(0).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(0)  << endl;
	cout << setw(10) << system.GetSpecies(1).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(1)  << endl;
	cout << setw(10) << system.GetSpecies(2).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(2)  << endl;
	cout << setw(10) << system.GetSpecies(3).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(3)  << endl;
	cout << setw(10) << system.GetSpecies(4).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(4)  << endl;
	cout << setw(10) << system.GetSpecies(5).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(5)  << endl;
	cout << setw(10) << system.GetSpecies(6).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(6)  << endl;
	cout << setw(10) << system.GetSpecies(7).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(7)  << endl;
	cout << setw(10) << system.GetSpecies(8).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(8)  << endl;
	cout << setw(10) << system.GetSpecies(9).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(9)  << endl;
	cout << setw(10) << system.GetSpecies(10).GetName() << ": " << setw(15) << system.GetMolesOfSpecies(10) << endl;
	cout << setw(10) << system.GetSpecies(11).GetName() << ": " << setw(15) << system.GetMolesOfSpecies(11) << endl;
	
	cout << endl;
	
	cout << log10(system.GetEquilibriumReactionQuotient(0)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(1)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(2)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(3)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(4)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(5)) << endl;
}

void TestBatchReactorWaterNaClCaCO3CO2()
{
	AqueousPhase aqueousPhase;
	
	aqueousPhase.AddSpecies("H2O(l)");
	aqueousPhase.AddSpecies("H+");
	aqueousPhase.AddSpecies("OH-");
	aqueousPhase.AddSpecies("Na+");
	aqueousPhase.AddSpecies("Cl-");
	aqueousPhase.AddSpecies("Ca++");
	aqueousPhase.AddSpecies("HCO3-");
	aqueousPhase.AddSpecies("NaCl(aq)");
	aqueousPhase.AddSpecies("CaCO3(aq)");
	aqueousPhase.AddSpecies("CO2(aq)");

	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
	aqueousPhase.SetActivityCoefficientModelSetschenow("CaCO3(aq)", 0.10);
	aqueousPhase.SetActivityCoefficientModelSetschenow("CO2(aq)", 0.20);
	
	GaseousPhase gaseousPhase;
	
	gaseousPhase.AddSpecies("H2O(g)");
	gaseousPhase.AddSpecies("CO2(g)");
	
	MineralPhase calcitePhase;
	
	calcitePhase.AddSpecies("Calcite");
	
	EquilibriumReaction eqReactionOH;
	
	eqReactionOH.SetNameOfReactant("OH-");
	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
	eqReactionOH.SetEquilibriumConstant({13.9951}, {25.0});
	
	EquilibriumReaction eqReactionNaClaq;
	
	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionNaClaq.SetEquilibriumConstant({0.7769}, {25.0});
	
	EquilibriumReaction eqReactionCaCO3aq;
	
	eqReactionCaCO3aq.SetNameOfReactant("CaCO3(aq)");
	eqReactionCaCO3aq.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
	eqReactionCaCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	eqReactionCaCO3aq.SetEquilibriumConstant({7.09941}, {25.0});
	
	EquilibriumReaction eqReactionCO2aq;
	
	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2aq.SetEquilibriumConstant({-6.35332}, {25.0});
	
	EquilibriumReaction eqReactionH2Og;
	
	eqReactionH2Og.SetNameOfReactant("H2O(g)");
	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
	eqReactionH2Og.SetEquilibriumConstant({1.49929}, {25.0});
	
	EquilibriumReaction eqReactionCO2g;
	
	eqReactionCO2g.SetNameOfReactant("CO2(g)");
	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2g.SetEquilibriumConstant({-7.82545}, {25.0});
	
	MineralReaction minReactionCalcite;
	
	minReactionCalcite.SetNameOfReactant("Calcite");
	minReactionCalcite.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
	minReactionCalcite.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	minReactionCalcite.SetEquilibriumConstant({1.84377}, {25.0});
	minReactionCalcite.SetSpecificSurfaceArea(9.8 * 100.09); // 9.8 cm2/g * 100.09 g/mol
	
	MineralReactionMechanism calciteMechanism;
	
	calciteMechanism.k  = pow(10.0, -3.48);
	calciteMechanism.Ea = 35.4;
	
	minReactionCalcite.AddMechanisms(calciteMechanism);

	System system;

	system.SetTemperature(25.0);
	system.SetPressure(100.0);
	
	system.AddPhase(aqueousPhase);
	system.AddPhase(gaseousPhase);
	system.AddPhase(calcitePhase);

	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "Ca++", "HCO3-"};
	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "CaCO3(aq)", "CO2(aq)", "H2O(g)", "CO2(g)"};
	vector<string> kineticSpecies = {"Calcite"};
	
	system.SetPrimarySpecies(primarySpecies);
	system.SetEquilibriumSpecies(equilibriumSpecies);
	system.SetKineticSpecies(kineticSpecies);
	
	system.AddEquilibriumReaction(eqReactionOH);
	system.AddEquilibriumReaction(eqReactionNaClaq);
	system.AddEquilibriumReaction(eqReactionCaCO3aq);
	system.AddEquilibriumReaction(eqReactionCO2aq);
	system.AddEquilibriumReaction(eqReactionH2Og);
	system.AddEquilibriumReaction(eqReactionCO2g);
	
	system.AddKineticReaction(minReactionCalcite);
	
	double C_H2O     = 55.508;	
	double C_H       = 1.0E-4;
	double C_Na      = 1.0E-4;
	double C_Cl      = 1.0E-4;
	double C_Ca      = 1.0E-4;
	double C_HCO3    = 1.0E-4;
	
	double C_OH      = 1.0E-4;
	double C_NaClaq  = 1.0E-4;
	double C_CaCO3aq = 1.0E-4;
	double C_CO2aq   = 1.0E-4;
	double C_H2Og    = 1.0E-4;
	double C_CO2g    = 1.0E-4;
	
	double C_Calcite  = 2.7;
	
	system.SetMolesOfSpecies({C_H2O, C_H, C_OH, C_Na, C_Cl, C_Ca, C_HCO3, C_NaClaq, C_CaCO3aq, C_CO2aq, C_H2Og, C_CO2g, C_Calcite});
	
	Speciation speciation;
	
	speciation.SetSystem(system);
	
	speciation.ImposeChargeBalance();
	speciation.ImposeCj(0, 55.508);	// H2O
	speciation.ImposeCj(2, 1);		// Na+
	speciation.ImposeCj(3, 1);		// Cl-
	speciation.ImposeCj(4, 1);		// Ca++
	speciation.ImposeCe(5, 2.0);	// CO2(g) (1 mole = 44 grams)
	
	speciation.Solve(1.0E-6, 100);
	
	cout << endl;
	
	cout << setw(10) << system.GetSpecies(0).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(0)  << endl;
	cout << setw(10) << system.GetSpecies(1).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(1)  << endl;
	cout << setw(10) << system.GetSpecies(2).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(2)  << endl;
	cout << setw(10) << system.GetSpecies(3).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(3)  << endl;
	cout << setw(10) << system.GetSpecies(4).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(4)  << endl;
	cout << setw(10) << system.GetSpecies(5).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(5)  << endl;
	cout << setw(10) << system.GetSpecies(6).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(6)  << endl;
	cout << setw(10) << system.GetSpecies(7).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(7)  << endl;
	cout << setw(10) << system.GetSpecies(8).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(8)  << endl;
	cout << setw(10) << system.GetSpecies(9).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(9)  << endl;
	cout << setw(10) << system.GetSpecies(10).GetName() << ": " << setw(15) << system.GetMolesOfSpecies(10) << endl;
	cout << setw(10) << system.GetSpecies(11).GetName() << ": " << setw(15) << system.GetMolesOfSpecies(11) << endl;
	
	cout << endl;
	
	cout << log10(system.GetEquilibriumReactionQuotient(0)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(1)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(2)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(3)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(4)) << endl;
	cout << log10(system.GetEquilibriumReactionQuotient(5)) << endl;
	
	BatchReactor batchReactor;
	
	batchReactor.SetSystem(&system);
	
	batchReactor.Solve(5, 1.0E-6, 1.0E-3, "BacthReactor-CO2.dat");
}

void TestBatchReactorWaterNaClCO2CalciteMagnesite()
{
	AqueousPhase aqueousPhase;
	aqueousPhase.AddSpecies("H2O(l)");
	aqueousPhase.AddSpecies("H+");
	aqueousPhase.AddSpecies("OH-");
	aqueousPhase.AddSpecies("Na+");
	aqueousPhase.AddSpecies("Cl-");
	aqueousPhase.AddSpecies("Ca++");
	aqueousPhase.AddSpecies("Mg++");
	aqueousPhase.AddSpecies("HCO3-");
	aqueousPhase.AddSpecies("NaCl(aq)");
	aqueousPhase.AddSpecies("HCl(aq)");
	aqueousPhase.AddSpecies("NaOH(aq)");
	aqueousPhase.AddSpecies("CaCO3(aq)");
	aqueousPhase.AddSpecies("MgCO3(aq)");
	aqueousPhase.AddSpecies("NaHCO3(aq)");
	aqueousPhase.AddSpecies("CO2(aq)");
	aqueousPhase.AddSpecies("CO3--");

	aqueousPhase.SetActivityCoefficientModelSetschenow("NaCl(aq)", 0.25);
	
	GaseousPhase gaseousPhase;
	gaseousPhase.AddSpecies("H2O(g)");
	gaseousPhase.AddSpecies("CO2(g)");
	
	MineralPhase calcitePhase;
	calcitePhase.AddSpecies("Calcite");
	
	MineralPhase magnesitePhase;
	magnesitePhase.AddSpecies("Magnesite");
	
	EquilibriumReaction eqReactionOH;
	
	eqReactionOH.SetNameOfReactant("OH-");
	eqReactionOH.SetNameOfProducts(vector<string>({"H2O(l)", "H+"}));
	eqReactionOH.SetStoichiometryOfProducts({1.0, -1.0});
	eqReactionOH.SetEquilibriumConstant({14.0013, 12.03594}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionNaClaq;
	
	eqReactionNaClaq.SetNameOfReactant("NaCl(aq)");
	eqReactionNaClaq.SetNameOfProducts(vector<string>({"Na+", "Cl-"}));
	eqReactionNaClaq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionNaClaq.SetEquilibriumConstant({0.7769, 0.473}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionHClaq;
	
	eqReactionHClaq.SetNameOfReactant("HCl(aq)");
	eqReactionHClaq.SetNameOfProducts(vector<string>({"H+", "Cl-"}));
	eqReactionHClaq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionHClaq.SetEquilibriumConstant({0.69754, 0.62}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionNaOHaq;
	
	eqReactionNaOHaq.SetNameOfReactant("NaOH(aq)");
	eqReactionNaOHaq.SetNameOfProducts(vector<string>({"Na+", "H2O(l)", "H+"}));
	eqReactionNaOHaq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	eqReactionNaOHaq.SetEquilibriumConstant({14.18, 12.19733}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionCaCO3aq;
	
	eqReactionCaCO3aq.SetNameOfReactant("CaCO3(aq)");
	eqReactionCaCO3aq.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
	eqReactionCaCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	eqReactionCaCO3aq.SetEquilibriumConstant({7.09941, 6.06247}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionMgCO3aq;
	
	eqReactionMgCO3aq.SetNameOfReactant("MgCO3(aq)");
	eqReactionMgCO3aq.SetNameOfProducts(vector<string>({"Mg++", "HCO3-", "H+"}));
	eqReactionMgCO3aq.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	eqReactionMgCO3aq.SetEquilibriumConstant({7.34394, 6.42957}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionNaHCO3aq;
	
	eqReactionNaHCO3aq.SetNameOfReactant("NaHCO3(aq)");
	eqReactionNaHCO3aq.SetNameOfProducts(vector<string>({"Na+", "HCO3-"}));
	eqReactionNaHCO3aq.SetStoichiometryOfProducts({1.0, 1.0});
	eqReactionNaHCO3aq.SetEquilibriumConstant({0.24474, 0.72845}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionCO2aq;
	
	eqReactionCO2aq.SetNameOfReactant("CO2(aq)");
	eqReactionCO2aq.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2aq.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2aq.SetEquilibriumConstant({-6.35332, -6.03078}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionCO3;
	
	eqReactionCO3.SetNameOfReactant("CO3--");
	eqReactionCO3.SetNameOfProducts(vector<string>({"HCO3-", "H+"}));
	eqReactionCO3.SetStoichiometryOfProducts({1.0, -1.0});
	eqReactionCO3.SetEquilibriumConstant({10.32666, 9.80904}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionH2Og;
	
	eqReactionH2Og.SetNameOfReactant("H2O(g)");
	eqReactionH2Og.SetNameOfProducts(vector<string>({"H2O(l)"}));
	eqReactionH2Og.SetStoichiometryOfProducts({1.0});
	eqReactionH2Og.SetEquilibriumConstant({1.49929, -0.05018}, {25.0, 100.0});
	
	EquilibriumReaction eqReactionCO2g;
	
	eqReactionCO2g.SetNameOfReactant("CO2(g)");
	eqReactionCO2g.SetNameOfProducts(vector<string>({"H2O(l)", "H+", "HCO3-"}));
	eqReactionCO2g.SetStoichiometryOfProducts({-1.0, 1.0, 1.0});
	eqReactionCO2g.SetEquilibriumConstant({-7.82545, -8.19835}, {25.0, 100.0});
	
	MineralReaction minReactionCalcite;
	
	minReactionCalcite.SetNameOfReactant("Calcite");
	minReactionCalcite.SetNameOfProducts(vector<string>({"Ca++", "HCO3-", "H+"}));
	minReactionCalcite.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	minReactionCalcite.SetEquilibriumConstant({1.84377, 0.99075}, {25.0, 100.0});
	minReactionCalcite.SetSpecificSurfaceArea(9.8 * 100.09); // 9.8 cm2/g * 100.09 g/mol
	
	MineralReactionMechanism calciteMechanism;
	
	calciteMechanism.k  = pow(10.0, -5.81);
	calciteMechanism.Ea = 23.5;
	
	minReactionCalcite.AddMechanisms(calciteMechanism);
	
	MineralReaction minReactionMagnesite;
	
	minReactionMagnesite.SetNameOfReactant("Magnesite");
	minReactionMagnesite.SetNameOfProducts(vector<string>({"Mg++", "HCO3-", "H+"}));
	minReactionMagnesite.SetStoichiometryOfProducts({1.0, 1.0, -1.0});
	minReactionMagnesite.SetEquilibriumConstant({2.03574, 0.4324}, {25.0, 100.0});
	minReactionMagnesite.SetSpecificSurfaceArea(9.8 * 84.31); // 9.8 cm2/g * 84.31 g/mol
	
	MineralReactionMechanism magnesiteMechanism;
	
	magnesiteMechanism.k  = pow(10.0, -9.34);
	magnesiteMechanism.Ea = 23.5;
	
	minReactionMagnesite.AddMechanisms(magnesiteMechanism);
	
	System system;

	system.SetTemperature(100.0);
	system.SetPressure(100.0);
	
	system.AddPhase(aqueousPhase);
	system.AddPhase(gaseousPhase);
	system.AddPhase(calcitePhase);
	system.AddPhase(magnesitePhase);

	vector<string> primarySpecies = {"H2O(l)", "H+", "Na+", "Cl-", "Ca++", "Mg++", "HCO3-"};
	vector<string> equilibriumSpecies = {"OH-", "NaCl(aq)", "HCl(aq)", "NaOH(aq)", "CaCO3(aq)", "MgCO3(aq)", "NaHCO3(aq)", "CO2(aq)", "CO3--", "H2O(g)", "CO2(g)"};
	vector<string> kineticSpecies = {"Calcite", "Magnesite"};
	
	system.SetPrimarySpecies(primarySpecies);
	system.SetEquilibriumSpecies(equilibriumSpecies);
	system.SetKineticSpecies(kineticSpecies);
	
	system.AddEquilibriumReaction(eqReactionOH);
	system.AddEquilibriumReaction(eqReactionNaClaq);
	system.AddEquilibriumReaction(eqReactionHClaq);
	system.AddEquilibriumReaction(eqReactionNaOHaq);
	system.AddEquilibriumReaction(eqReactionCaCO3aq);
	system.AddEquilibriumReaction(eqReactionMgCO3aq);
	system.AddEquilibriumReaction(eqReactionNaHCO3aq);
	system.AddEquilibriumReaction(eqReactionCO2aq);
	system.AddEquilibriumReaction(eqReactionCO3);
	system.AddEquilibriumReaction(eqReactionH2Og);
	system.AddEquilibriumReaction(eqReactionCO2g);
	
	system.AddKineticReaction(minReactionCalcite);
	system.AddKineticReaction(minReactionMagnesite);
	
	vector<double> C(system.GetNumOfSpecies(), 1.0E-7);
	
	C[system.GetIndexOfSpecies("H2O(l)")]    = 55.508;
	C[system.GetIndexOfSpecies("Na+")]       = 1.0;
	C[system.GetIndexOfSpecies("Cl-")]       = 1.0;
	C[system.GetIndexOfSpecies("Ca++")]      = 1.0E-7;
	C[system.GetIndexOfSpecies("Mg++")]      = 1.0E-7;
	C[system.GetIndexOfSpecies("HCO3-")]     = 2.0;
	C[system.GetIndexOfSpecies("CO2(g)")]    = 2.0; // Moleweight: 44.01 g/mol
	C[system.GetIndexOfSpecies("Calcite")]   = 2.0; // Moleweight: 100.09 g/mol, volumic mass: 2709.89 kg/m3
	C[system.GetIndexOfSpecies("Magnesite")] = 2.0; // Moleweight: 84.31 g/mol, volumic mass: 3009.28 kg/m3
	
	system.SetMolesOfSpecies(C);
	
	Speciation speciation;
	
	speciation.SetSystem(system);
	
	speciation.ImposeChargeBalance();
	speciation.ImposeCj(0, 55.508);	// H2O
	speciation.ImposeCj(2, 1.0);	// Na+
	speciation.ImposeCj(3, 1.0);	// Cl-
	speciation.ImposeCj(4, 1.0E-7);	// Ca++
	speciation.ImposeCj(5, 1.0E-7);	// Mg++
	speciation.ImposeCe(10, 2.0);	// CO2(g) (1 mole = 44 grams)
	
	speciation.Solve(1.0E-6, 100);
	
	cout << endl;
	
	for(unsigned iSpecies = 0; iSpecies < system.GetNumOfSpecies(); ++iSpecies)
		cout << setw(10) << system.GetSpecies(iSpecies).GetName()  << ": " << setw(15) << system.GetMolesOfSpecies(iSpecies)  << endl;
	
	cout << endl;
	
	for(unsigned iEquilibrium = 0; iEquilibrium < system.GetNumOfEquilibriumSpecies(); ++iEquilibrium)
		cout << "K[" << system.GetSpecies(system.GetIndexOfEquilibriumSpecies(iEquilibrium)).GetName()  << "]: " << setw(15) << log10(system.GetEquilibriumReactionQuotient(iEquilibrium)) << endl;
	
	BatchReactor batchReactor;
	
	batchReactor.SetSystem(&system);
	
	double hour = 3600.0;
	double day  = 24 * hour;
	
	batchReactor.Solve(5 * hour, 0.1, 0.1, "BacthReactor.Calcite-Magnesite-CO2-Water-NaCl.dat");
}

int main()
{
	//TestActivityCoefficients_NaCl();
	
	//TestSpeciationWaterNaCl();
	
	//TestSpeciationWaterNaClCaCO3();
	
	//TestBatchReactorWaterNaClCaCO3();
	
	//TestBatchReactorWaterNaClCaCO3CO2();
	
	//TestBatchReactorWaterNaClCO2CalciteMagnesite();
	
	//TestSpeciationWaterCO2();
	
	TestSpeciationWaterNaClCO2();
	
	//TestSpeciationWaterNaClCaCO3CO2();
}

